Tkachenko modes in a superfluid Fermi gas at unitarity 
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We calculate the frequencies of the Tkachenko oscillations of a vortex lattice in a harmonically 
trapped superfluid Fermi gas. We use the elasto-hydrodynamic theory by properly accounting for 
the elastic constants, the Thomas-Fermi density profile of the atomic cloud, and the boundary 
conditions. Thanks to the Fermi pressure, which is responsible for larger cloud radii with respect 
to the case of dilute Bose-Einstein condensed gases, large vortex lattices are achievable in the 
unitary limit of infinite scattering length, even at relatively small angular velocities. This opens 
the possibility of experimentally observing vortex oscillations in the regime where the dispersion 
relation approaches the Tkachenko law for incompressible fluids and the mode frequency is almost 
comparable to the trapping frequencies. 
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Two component Fermi gases are very versatile sys- 
tems, offering the possibility of exploring both fermionic 
and bosonic superfluidity. By exploiting Feshbach res- 
onances to tune the interaction strength between dif- 
ferent spin species, it is possible to realize ultracold 
gases of fermionic atoms and Bose-Einstein condensates 
(BECs) of dimers as well (for a recent review see Ref. [l|). 
Moreover, at resonance, where the interspecies scattering 
length diverges, fermionic atoms enter a new strongly cor- 
related regime, the so-called unitary limit, in which the 
standard Bardeen-Cooper-Schrieffer (BCS) theory is not 
applicable. Here fermionic superfluidity is significantly 
enhanced, due to the stronger role played by the inter- 
actions. Vortex lattices have been already observed in 
this novel phase of matter Q. The observation of vor- 
tices in Fermi gases represents a particularly important 
achievement, since the unambiguous detection of super- 
fluidity in these systems is less straightforward than in 
the corresponding case of Bose-Einstein condensates. 

In this paper we study the Tkachenko modes of the vor- 
tex lattice in a harmonically trapped Fermi gas at unitar- 
ity. These oscillations, originally studied in Ref. [3j for 
incompressible superfluids, correspond to shear distor- 
tions of the lattice planes and carry much lower energies 
than usual hydrodynamic modes. The investigation of 
Tkachenko modes in Bose-Einstein condensates has been 
largely pursued on both the theoretical 0, @, H, 0, H, 
and experimental [l(| [H| sides. An appealing reason to 
study vortex modes in Fermi gases at unitarity is given by 
the larger number of vortices N v achievable for a given 
rotation rate. As we shall discuss in the final part of 
this paper, unitary fcrmions of 6 Li yield a significant en- 
hancement (almost an order of magnitude) in the vortex 
number compared to bosonic atoms of 87 Rb. This is due 
to the large expansion of the cloud radius produced by 
the Fermi pressure. This effect provides the possibility 
of studying Tkachenko modes in a more systematic way 



than in the bosonic case, including an easier detection 
of modes with more than one radial node and the ob- 
servation of Tkachenko modes at relatively small angular 
velocity f2, where the incompressible regime of the dis- 
persion relation is exploited. 

A first insight into the problem can be gained from the 
frequency spectrum of the Tkachenko oscillations of an 
infinite vortex lattice fijl 
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where k is the wave vector, c s is the usual sound velocity, 
ct = yj kCI/8tt is the Tkachenko velocity, and the limit 
c s 3> ct is assumed. Here k = h/M is the quantum of cir- 
culation of a single vortex line, where M = m for bosons 
and M = 2m for fermions, m being the mass of the par- 
ticles. The incompressible limit of Eq. JT]) takes place for 
c s k ^> f2 and corresponds to the original Tkachenko dis- 
persion law lu — crk. In the opposite limit c s k <C SI one 
finds the quadratic dispersion relation u> — CTC s k 2 /2fi, 
characterizing the compressible limit of the spectrum. In 
a non-uniform configuration the actual values of k are 
inversely proportional to the radius of the cloud and are 
fixed by the proper procedures of discretization that will 
be discussed later. In the experiments with BECs it is 
difficult to realize large vortex lattices at small angular 
velocities and the condition c s k 3> SI is hence hard to ful- 
fill. Vice versa, we will show that in the case of fermions 
at unitarity the crossover from incompressible to com- 
pressible behavior can be more easily investigated. 

In our calculation, we adapt the two-dimensional (2D) 
elasto-hydrodynamic treatment developed by Sonin [7| to 
the fermionic case. We consider two distinct density pro- 
files in the Thomas-Fermi (TF) approximation: (i) the 
one corresponding to the cylinder geometry as consid- 
ered in Ref. [7| and (ii) the 2D column density obtained 
by integrating a 3D TF cloud along the axial direction 
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(hereafter also called pancake geometry). The cylinder 
geometry is an idealized configuration whose interest lies 
in the exact decoupling between the radial and axial dy- 
namics. However, experimental situations are expected 
to be better described by the pancake geometry, which 
takes into account the 3D features of the inhomogeneous 
profile. While in the full 3D geometry the decoupling of 
dynamics is not exact, the axial motion is expected to 
be negligible for Tkachenko modes as long as vortices are 
straight. The use of the column density is hence expected 
to be a valuable approximation for pancake shape clouds 
which do not exhibit significant vortex bending. 

One of the major differences between Bose and Fermi 
gases lies in the equation of state, being dominated in the 
Fermi case by the quantum pressure effect. In the TF 
approximation, the equation of state can be expressed as 
u cx n 1 in both cases, but with different values of the 
polytropic index 7. Here u is the chemical potential and 
n is the number density of particles. In 3D, 7 = 1 for 
bosons and 7 = 2/3 for fermions at unitarity pj|], yield- 
ing a different TF density profile in the two cases. For a 
vortex lattice rotating at an angular velocity f2 in a har- 
monic trap with radial and axial frequencies lo± and u> z , 
the coarse-grained equilibrium density profile is given by 
n(r) cx [/Jo — y(r)] 1 / 7 , where fig is the chemical poten- 
tial for the trapped case, V(r) = rn{G)\r\ + uj 2 z 2 )/2 is 
the effective potential taking into account the centrifugal 
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For a TF configuration with pancake geometry we may 
reduce the problem to two dimensions by integrating out 
the z coordinate; i.e., we employ the column density 
n2T>(r±) = J dzn(r) as an effective 2D density profile. 
One hence finds 712d(»"_l) oc (fio — mtZ;^r^/2) 1 / 7off with 
7 c ff = 27/(2 + 7) corresponding to the equation of state 
a cx n 7off in the uniform 2D configuration. Note that the 
effective 2D polytropic index 7 e ff always differs from the 
one calculated in the cylindrical geometry, where lo z = 
and 7eff = 7- It is also worth pointing out that j e g has 
the same value for fermionic systems in the cylindrical 
geometry and for bosonic systems in the pancake geom- 
etry (see Table Uj. In the remainder of this paper, we 
treat the problem only in two dimensions and drop the 
subscript "_L" of r± and the subscript "2D" of ri2D- 

Macroscopic manifestations of superfluidity at zero 
temperature can be studied by resorting to the super- 
fluid hydrodynamic equations. As long as we study 
the dynamics of vortices on length scales much larger 
than the intervortex distance, we can employ the elasto- 
hydrodynamic theory [l^, Q EH • Here the microscopic 
density and velocity fields are substituted by averaged 
quantities and the restoring force of the lattice is taken 
into account by adding an elastic energy term to the usual 
superfluid hydrodynamic energy functional [lr| . 

The elastic energy for a triangular lattice is given by 
E c \ = J dr£ c \, where the energy density £ e \ in the rotat- 
ing frame is £ c \ = 2Ci(V ■ e) 2 + C 2 [(de x / dx — de y /dy) 2 + 
(de x /dy + de y /dx) 2 ]. Here e is the vortex displacement 



TABLE I: Effective polytropic index 7 e ff in 2D for the bosonic 
and fermionic cases. For the cylindrical geometry, y e g — 7, 
and for the pancake geometry, y e ff — 27/(2 + 7). The values 
of the parameter a characterizing the dispersion law in the 
incompressible limit (see text) are also reported. 



geometry 


bosons 


fermions 


cylinder 


1 (a = 5.43) 


2/3 (a = 5.59) 


pancake 


2/3 (a = 5.59) 


1/2 (a = 5.75) 



field, and the coefficients C\ and C2 are the compres- 
sional and shear modulus respectively [Til ] , corresponding 
to second derivatives of the energy density with respect 
to lattice distortions. They have to be calculated from 
the (microscopic) energy functional evaluated in the ro- 
tating frame [151 ]. 

In the Fermi case at unitarity the core size £ of vortices 
is of the order of the intcrparticle distance. Consequently, 
unless one works extremely close to the centrifugal limit 
[l||, £ is much smaller than the intervortex distance, 
which is of the order of the Wigne r-Seitz radius of the 
vortex lattice cell, namely, In = \J k/2itVL = ^h/MVl. 
This relation corresponds to the usual vortex density 
n v = 1/7tZj 2 = 2fl/ k obtained from the condition of quan- 
tized circulation. The elastic coefficients can then be cal- 
culated in the small core limit £ <C In, equivalent to the 
Thomas- Fermi condition ^0 3> Ml [20J. In this regime, 
the well-established result C2 = —C\ = nmntt / '16tt holds 
fl2L IT5L l2ll] . Hence, for bosons C2 = — G\ = nhfl/8, while 
for fermions Ci — —C\ — nHCt/16, i.e., in the latter case 
n is replaced by the density of pairs, n/2. This result can 
be understood by observing that fermionic pairs play the 
same role as bosonic molecules from the point of view of 
superfluidity. 

We are finally ready to write the linearized elasto- 
hydrodynamic equations. In the rotating frame they take 
the form 



d_ 

di 



Sn - 



— dv + 2nAdv + V 

at to 



V- (n 5v) 
Sfj, F el 



dt 



5v + 2n Ae + V 



777770 
5 ft 



(2) 
(3) 
(4) 



where e = de/dt and the variation of the local chemical 
potential, 5fi, can be expressed in terms of the local sound 
velocity as 5/i/m — (dfi/dn)5n/m = (c 2 /no)Sn. By com- 
bining Eqs. §5§ and (H]) one also finds 2mnoflA(Sv — e) = 
F c \. Here 77 = 77 (r) is the equilibrium density, while 
Sn and Sv are the density and velocity perturbations. 
The elastic force is given by (F e i)j — d<Tjk/dxk, where 
<jjk = SE C \/Sujk = d£ c \/dujk is the stress tensor defined 
in terms of the strain tensor ?jjfc = (dej / dxk+dek/dxj) /2 
Note that, while in the case of BECs for the cylin- 
drical geometry the ratio c 2 s /no is constant and commutes 
with the gradient in Eqs. © and (]3|), this is not the case 
for 7eg 7^ 1- 
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FIG. 1: (Color online) Plot of the frequency w of the first 
and the second axisymmetric Tkachenko modes as a function 
of the angular velocity. The green dashed line corresponds to 
pancake bosons ( 7e ff — 2/3), while the red solid line to pan- 
cake fermions (y c f{ = 1/2). For each mode the incompressible 
limit Wi n o = iooa is shown for comparison. 



We are interested in the axisymmetric Tkachenko 
modes. In polar coordinates all the physical quantities 
are then independent of the azimuthal angle <p and we 
can reduce to ordinary differential equations with respect 
to r. For large vortex numbers one can approximate 
£0 ~ Svcf, (see Ref. [7]) and the above linearized equa- 
tions take the simplified form 
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where Sv r and 5v^ are the radial and azimuthal compo- 
nents of Sv, and ct = ^nfl/Sir = Iq,£1/2. 

We solve Eqs. ([5} and {B| with proper boundary con- 
ditions at the cloud center and at the cloud radius 
R = y / 2fj 1 /muj 2 _ using the shooting method. At r = 
the velocity must vanish: Sv r (0)_= 5v<p(0) = 0. Follow- 
ing the same arguments of Ref. [7[ , we obtain at r = R 
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In Fig. [T] we present our predictions for the frequen- 
cies of the two lowest Tkachenko modes at unitarity. In 
the numerical calculation it is natural to solve Eqs. ([5]) 
and ^ for the dimensionless ratio oj/ujq, where luq = 
ct/R. Since the discretized values of k are proportional 
to l/R the frequency luq is basically the analogue of 
the incompressible limit crk of Eq. (fTJ). Comparison 
with the homogeneous spectrum of Eq. (fT]) then shows 
that uj/ujo is the analogue of c s k/ 'y '4ft 2 + c 2 k 2 , i.e., of 
the factor responsible for compressibility effects. The 
dependence of lu on the angular velocity can be ex- 
pressed analytically. The radius of the cloud is given 



by R = R /(l - 
rotating cloud, f2 = 



v , where Rq is the radius of the non- 
(l/ux, and v = [2( 7cf f + 1)]" 1 . Then, 
- tt 2 Y, where A Q = v / m/4Ma ± /R Q 

and a± = h/mu>±. The full dependence of u> on is 
hence displayed by plotting the quantity u>/Aquj± which 
does not depend any longer on Rq, but only on the ratio 
n/w±. In the figure we also report the predictions for a 
dilute Bose-Einstein condensed gas in the same pancake 
configuration. The differences are actually very small, 
showing that, when expressed in the units of Fig. [T] the 
results for the Tkachenko frequencies do not depend in 
an appreciable way on the actual form of the equation of 
state which instead can significantly affect the value of 
the TF radius and hence of A . 

One can get a qualitative understanding of the depen- 
dence of u! on Q in the trapped case looking at the homo- 
geneous dispersion relation (fT]). To this purpose, one has 
to evaluate c s and k in the finite size system. A simple 
estimate for the sound velocity is given by its value at 
the center of the cloud, i.e., c 2 ~ 7 c ff/i/?7i = r y c ff£j 2 _R 2 /2. 
The effective wave vector can instead be quantized pro- 
portionally to l/R, where the proportionality coefficient 
can be extracted by comparing Eq. {1} with the numer- 
ical calculation. For a general value of fi, the latter re- 
sult would depend on the chosen estimate for c s . How- 
ever, in the fl — > limit Eq. (fT]) becomes independent 
of c s and one can unambiguously define k = a/R with 
a = limr2_*o w/cjq. The quantized value of k can then 
be used to estimate the full O dependence using Eq.([T]). 
This procedure was first exploited by Baym [5( in the 
case of BECs and yields a rather accurate estimate of 
the whole dispersion law [22j]. We find that the value of 
a exhibits only a weak dependence on 7 and on the ge- 
ometry employed (see Table [J. For the lowest Tkachenko 
mode we find the value a — 5.75 for fermions in the pan- 
cake geometry to be compared with the result a = 5.59 
holding for bosons in the same geometry (or for fermions 
in the cylindrical geometry) [23|]. Note also that the in- 
compressible limit Wi nc = u>qcx of the dispersion relation 
reproduces the full dispersion with high accuracy up to 
values (l — 0.4 ~ 0.5 of the angular velocity (see Fig. [1]). 

Let us now discuss the dependence of the number of 
vortices N v on the experimental parameters. This is 
important in order to determine the physical regimes 
of £1/uj± achievable in practice. The number of vor- 
tices in the TF approximation is N v = n v irR 2 . For 
a harmonically trapped Fermi gas at unitarity one has 
H R = a ± [2VTT^(3A^A) 1 / 3 ] 1 / 2 (l - ft 2 )- 1 / 3 , where 
/3 ~ —0.6 is a universal dimensionless parameter account- 
ing for the role of interactions [25[, N is the number 
of particles, and A = uj z /uj± is the trap aspect ratio. 
Thus, N v (fermions) = iyTTpflil - n 2 )- 2 / 3 (3N\)^ 3 . 
Using the corresponding expression for a Bose-Einstein 
condensed gas we find 



N v (fermions) 4yi + (3 
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where the subscripts refer to bosons and fermions and 
a is the s-wave scattering length. For (NX) / = 
(NX)b the above ratio becomes 4-yT + (3 [3iVA(l — 
£l 2 )]~ 1 / 15 (5a/a±) b 2 ^ 5 and the small exponent —1/15 im- 
plies that the result is practically insensitive to the value 
of fl, N, and A. On the other hand, in typical BECs one 
has a/a± < 1, so that N v (fermions) /N v (bosons) can be 
significantly larger than 1 [26| . 

In the JILA experiments on BECs one has A = 0.63, 
o~5nm, = 3.74 fj,m, and N ~ 10 6 \u%. This yields 
N v ( 87 Rb) = 43.7fi(l - ft 2 )- 3 / 5 . In the MIT experi- 
ments, at unitarity one has A = 0.40 and N ~ 10 6 Q. 
Then AT„( 6 Li)//V„( 87 Rb) ~ 6.1(1 - ft 2 )" 1 / 15 , which cor- 
responds to the significant gain in the number of vortices 
mentioned in the introduction (27j . 

The increase of N v at unitarity allows to measure 
Tkachenko modes at relatively small angular velocities. 
For example, a fermionic cloud with the above parame- 
ters easily contains more than 100 vortices at O = 0.4, 
deeply in the incompressible region of the spectrum. The 



maximum of the frequency of the lowest Tkachenko mode 
takes place at ft ~ 0.55 and corresponds to ui/u>± ~ 0.13 
for the same parameters. This value is quite larger 
than the highest frequency observed in 87 Rb experiments, 
namely, uj/oj^ = 0.023 at O = 0.84 [lo|. This pro- 
vides promising pers pec tives for precision measurements 
of Tkachenko modes 28] , for a direct determination of the 
quantum of circulation in Fermi superfluids, and for the 
general investigation of the incompressible regime of the 
dispersion relation, whose quantitative analysis in experi- 
ments is a long standing question in the field of superfluid 
systems. 
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